SPARSITY  ADAPTIVE  MATCHING  PURSUIT  ALGORITHM  FOR  PRACTICAL 

COMPRESSED  SENSING 


Thong  T.  Do1,  Lu  GarT,  Nam  Nguyen t  and  True  D.  Tran  '  * 


t  Department  of  Electrical  and  Computer  Engineering 
The  Johns  Hopkins  University 
* School  of  Engineering  and  Design 
Brunei  University,  UK 


ABSTRACT 

This  paper  presents  a  novel  iterative  greedy  reconstruction  algorithm 
for  practical  compressed  sensing  (CS),  called  the  sparsity  adaptive 
matching  pursuit  (SAMP).  Compared  with  other  state-of-the-art 
greedy  algorithms,  the  most  innovative  feature  of  the  SAMP  is  its 
capability  of  signal  reconstruction  without  prior  information  of  the 
sparsity.  This  makes  it  a  promising  candidate  for  many  practical 
applications  when  the  number  of  non-zero  (significant)  coefficients 
of  a  signal  is  not  available.  The  proposed  algorithm  adopts  a  sim¬ 
ilar  flavor  of  the  EM  algorithm,  which  alternatively  estimates  the 
sparsity  and  the  true  support  set  of  the  target  signals.  In  fact,  SAMP 
provides  a  generalized  greedy  reconstruction  framework  in  which 
the  orthogonal  matching  pursuit  and  the  subspace  pursuit  can  be 
viewed  as  its  special  cases.  Such  a  connection  also  gives  us  an  in¬ 
tuitive  justification  of  trade-offs  between  computational  complexity 
and  reconstruction  performance.  While  the  SAMP  offers  a  compara¬ 
bly  theoretical  guarantees  as  the  best  optimization-based  approach, 
simulation  results  show  that  it  outperforms  many  existing  iterative 
algorithms,  especially  for  compressible  signals. 

Index  Terms —  Sparsity  adaptive,  greedy  pursuit,  compressed 
sensing,  compressive  sampling,  sparse  reconstruction 

1.  INTRODUCTION 

Compressed  sensing  (CS)  [1]  has  gained  increased  interests  over  the 
past  few  years.  Suppose  that  x  is  a  length-lV  signal.  It  is  said  to 
be  Iv -sparse  (or  compressible)  if  x  can  be  well  approximated  using 
K  <C  N  coefficients  under  some  linear  transform.  According  to  the 
CS  theory,  such  a  signal  can  be  acquired  through  the  following  linear 
random  projections: 

y  =  $x  +  e,  (1) 

where  y  is  the  sampled  vector  with  M  <C  N  data  points,  $  rep¬ 
resents  an  M  x  N  random  projection  matrix  and  e  is  the  acquisi¬ 
tion  noise.  The  CS  framework  is  attractive  as  it  implies  that  x  can 
be  faithfully  recovered  from  only  M  =  0(K  log  N)  samples  [1], 
suggesting  the  potential  of  significant  cost  reduction  in  digital  data 
acquisition. 

Although  the  encoding  process  is  simply  linear  projection,  the 
reconstruction  requires  some  non-linear  algorithms  to  find  the  spars¬ 
est  signal  from  the  measurements.  One  challenging  question  in  the 
CS  research  is  the  development  of  fast  reconstruction  algorithm 
with  reliable  accuracy  and  (nearly)  optimal  theoretical  performance. 
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Among  existing  reconstruction  algorithms,  the  famous  basis  pur¬ 
suit  (BP)  [2]  aims  at  the  l\  minimization  using  linear  programming 
(LP).  While  it  requires  a  minimal  number  of  measurements,  its 
high  computational  complexity  may  prevent  it  from  practical  large- 
scale  applications.  Several  fast  convex  relaxation  algorithms  have 
been  proposed  to  solve  or  approximate  the  solution  of  BP,  e.g..  the 
gradient  projection  method  in  [3], 

Another  popular  class  of  sparse  recovery  algorithms  is  based  on 
the  idea  of  iterative  greedy  pursuit.  The  earliest  ones  include  the 
matching  pursuit  and  orthogonal  matching  pursuit  (OMP)  [4],  Their 
successors  include  the  stagewise  OMP  (StOMP)  [5]  and  the  regular¬ 
ized  OMP  (ROMP)  [6].  The  reconstruction  complexity  of  these  al¬ 
gorithms  is  around  0(KMN),  which  is  significantly  lower  than  the 
BP  methods.  However,  they  require  more  measurements  for  perfect 
reconstruction  and  they  lack  provable  reconstruction  quality.  More 
recently,  greedy  algorithms  such  as  the  subpace  pursuil(SP)  [7]  and 
the  compressive  sampling  matching  pursuit  (CoSaMP)  [8]  have  been 
proposed  by  incorporating  the  idea  of  backtracking.  They  offer  com¬ 
parable  theoretical  reconstruction  quality  as  that  of  the  LP  methods 
and  low  reconstruction  complexity.  However,  both  the  SP  and  the 
CoSAMP  assume  that  the  sparsity  K  is  known,  whereas  K  may  not 
be  available  in  many  practical  applications. 

In  this  paper,  we  propose  a  new  greedy  algorithm  called  the 
sparsity  adaptive  matching  pursuit  (SAMP)  for  blind  signal  recovery 
when  K  is  unknown.  SAMP  is  a  generalization  of  existing  greedy 
algorithms  as  both  the  OMP  and  the  SP  can  be  viewed  as  its  special 
cases.  It  follows  the  “divide  and  conquer”  principle  through  stage  by 
stage  estimation  of  the  sparsity  level  and  the  true  support  set  of  the 
target  signals.  The  SAMP  offers  a  comparably  theoretical  guaran¬ 
tees  as  the  best  optimization-based  approach.  Its  numerical  results 
are  even  more  attractive  as  it  outperforms  all  of  the  above-mentioned 
algorithms  in  extensive  simulations,  including  the  h  -minimization. 
The  rest  of  this  paper  is  organized  as  follows.  Section  2  depicts 
the  big  picture  of  above  mentioned  greedy  pursuit  algorithms  and 
presents  the  main  motivation  of  this  work.  While  detailed  descrip¬ 
tions  of  the  proposed  SAMP  algorithm  are  provided  in  Section  3, 
Section  4  presents  the  theoretical  analysis  of  exact  recovery  and  sta¬ 
bility.  Finally,  simulation  results  and  discussion  are  shown  in  Section 
5,  followed  by  the  conclusion  in  Section  6. 

2.  REVIEW 

This  section  presents  a  summary  of  existing  greedy  recovery  al¬ 
gorithms.  They  were  grouped  into  three  categories,  as  show  in 
Fig.  l(a)-(c).  Here,  in  the  fc-th  iteration,  r*  and  Fr  represent  the 
residue  and  the  estimated  signal’s  support  (called  finalist),  respec- 
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tively.  And  in  Fig.  1(c),  Ck  corresponds  to  the  candidate  set  of  the 
SP/CoSAMP  algorithms.  Among  these  algorithms,  the  OMP,  the 
STOMP  and  the  ROMP  adopt  a  bottom-up  approach  by  sequentially 
adding  the  support  set  of  x.  On  the  other  hand,  the  SP  and  the 
CoSaMP  use  a  top-down  approach  to  iteratively  refine,  rather  than 
expand  the  finalist. 

As  can  be  readily  seen,  while  the  OMP  and  the  StOMP  in 
Fig.  1(a)  use  only  one  test,  the  ROMP  in  Fig.  1(b)  uses  two  tests  to 
add  one  or  several  coordinates  to  the  finalist.  Specifically,  the  OMP 
adopts  the  maximal  correlation  test  and  in  each  iteration,  only  one 
candidate  is  added.  The  test  of  StOMP  follows  the  matched  filtering 
and  hard  thresholding  principles  to  choose  a  subset  of  coordinates  of 
atoms.  While  for  the  ROMP,  it  applies  a  preliminary  test  and  a  final 
test  to  build  the  finalist.  The  preliminary  test  is  quite  similar  to  that 
of  the  StOMP,  and  the  final  test  is  designed  to  keep  the  largest  subset 
of  those  coordinates  whose  atoms’  correlation  differ  in  magnitude 
by  at  most  a  factor  of  two.  For  these  bottom-up  methods,  the  finalist 
is  updated  at  the  end  of  iteration  by  union  of  the  new  discovered 
coordinates  and  the  finalist  in  the  previous  iteration.  Then,  the  ob¬ 
servation  residual  is  also  updated  by  subtracting  the  observation  data 
from  its  projection  onto  the  subspace  spanned  by  the  atoms  in  the 
finalist.  This  step  is  also  termed  as  orthogonalization  to  ensure  that 
the  observation  residual  is  always  orthogonal  to  atoms  in  the  finalist. 

Top-down  approaches  such  as  the  SP  and  the  CoSaMP  also  use 
two  different  tests  in  each  iteration.  But,  the  size  of  then  finalist  is 
kept  fixed  (and  equal  to  K  -  the  sparsity  of  input  signal).  In  partic¬ 
ular,  the  preliminary  test  is  quite  similar  to  that  of  the  ROMP  and 
the  final  test  is  designed  to  be  more  subtle  and  thus  more  reliable. 
After  the  preliminary  test,  a  candidate  list  is  created  by  union  of  the 
short  list  and  the  finalist  in  the  previous  iteration.  The  final  test  first 
solves  a  least  square  solution  and  then  choose  from  the  candidate  list 
a  subset  of  K  coordinates  that  are  corresponding  to  largest  entries  in 
magnitude  of  the  least  square  solution.  This  subset  of  coordinates 
serves  as  the  finalist.  The  observation  residual  is  finally  updated  in  a 
way  similar  to  that  of  above-mentioned  algorithms.  Compared  with 
the  bottom-up  greedy  algorithms,  the  remarkable  innovation  of  the 
SP  and  the  CoSaMP  is  the  backtracking  technique  in  their  final  test, 
which  enables  the  algorithms  to  remove  wrong  coordinates  added 
in  the  previous  iteration.  Among  existing  greedy  approaches,  only 
the  SP  and  the  CoSaMP  have  a  strong  theoretical  guarantee  com¬ 
parable  with  that  of  the  BP  (li-minimization).  Besides,  they  can 
operate  when  the  measurements  are  inaccurate  and/or  the  signal  is 
not  strictly  sparse. 

Despite  their  low  complexity,  all  greedy  pursuit  algorithms  in 
Fig.  1  require  the  sparsity  K  as  a  prior  for  exact  recovery.  How¬ 
ever,  in  practical  CS,  this  piece  of  information  is  often  not  avail¬ 
able.  For  example,  most  natural  image  signals  are  only  compressible 
(rather  than  strictly  sparse)  under  a  de-correlating  transform  (e.g., 
the  wavelet).  The  sparsity  K  (i.e.,  no.  of  significant  coefficients)  for 
these  signals  could  not  be  well-defined,  let  alone  be  known.  Some 
existing  algorithms  might  be  modified  to  handle  this  case.  For  ex¬ 
ample,  we  could  change  the  halting  condition  in  the  OMP,  i.e.,  iter¬ 
ating  until  the  energy  of  residual  is  smaller  than  a  certain  threshold. 
However,  it  is  not  known  if  this  modified  OMP  has  any  theoreti¬ 
cal  guarantee  of  exact  recovery  or  stability  yet.  In  another  way,  we 
might  want  to  guess  the  value  of  K  for  the  SP  or  the  CoSaMP.  How¬ 
ever,  it  would  either  eliminate  the  ability  of  exact  recovery  if  we 
underestimate  K  or  significantly  degrade  both  accuracy  and  robust¬ 
ness  of  the  algorithm  if  we  overestimate  it.  The  following  experi¬ 
ment  demonstrates  the  performance  degradation  of  the  SP  algorithm 
when  K  is  over-estimated.  Here,  x  is  a  uniformly  Gaussian  ran¬ 
dom  sparse  signals  with  length  of  N  =  256  and  the  sparsity  of 
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Fig.  1.  Conceptual  diagrams  of  (a)  the  OMP  and  the  STOMP;  (b) 
the  ROMP;  (c)  the  SP  and  the  CoSAMP. 


Fig.  2.  Probability  of  exact  recovery  vs.  estimation  of  sparsity  in  the 
SP  algorithm.  Here,  the  original  signal  has  length  of  N  =  256  with 
K  =  20  non-zero  entries. 

I\  =  20.  The  sensing  matrix  <I>  is  the  partial  Fourier  transform 
(PFFT)  [1],  500  simulations  were  conducted  for  each  pair  of  esti¬ 
mated  sparsity  K  =  (20, 30, 40,  50)  and  the  number  of  measure¬ 
ments  M  =  (50,  60,  70,  80, 90, 100).  Fig.  2  shows  the  curves  of  the 
probability  of  exact  recovery  vs.  estimation  of  the  sparsity.  One  can 
see  clearly  that  the  performance  of  the  SP  algorithm  drops  quickly  if 
the  estimated  sparsity  K  is  far  from  the  truth.  We  also  found  that  the 
CoSAMP  algorithm  shows  a  similar  performance  degradation.  Tak¬ 
ing  this  fact  into  account,  we  aim  to  develop  a  new  greedy  algorithm 
for  blind  recovery  when  the  sparsity  K  is  not  available. 

3.  SPARSITY  ADAPTIVE  MATCHING  PURSUIT 
3.1.  Algorithm  description 

Note  that  top-down  methods  such  as  the  SP  and  the  CoSAMP  are 
likely  to  identify  the  true  support  set  more  accurately  thanks  to  their 
backtracking  strategy.  On  the  other  hand,  bottom-up  approaches 
such  as  the  OMP  suggest  a  possible  solution  to  estimate  the  value 
of  K  by  moving  forward  step  by  step.  Following  these  observa- 
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Fig.  3.  A  Conceptual  Diagram  of  the  Proposed  SAMP 


Input:  Sampling  matrix  <f>.  Sampled  vector  y,  step  size  s; 
Output:  A  A'-sparse  approximation  x  of  the  input  signal; 

Initialization: 

x  =  0  {  Trivial  initialization} 

ro  =  y  {  Initial  residue} 

Fo  =  0  {  Empty  finalist} 

X  =  s  {  Size  of  the  finalist  in  the  first  stage} 
k  =  1  {  Iteration  index} 

j  —  1  {  Stage  index} 

repeat 

Sk  =  Max(\&*rk-i\,X )  {  Preliminary  Test} 

Ck  =  Fk- 1  U  Sk  {  Make  Candidate  List} 

F  =  Max(\&cky\,I)  {  Final  Test} 
r  =  y  —  {  Compute  Residue} 

if  halting  condition  true  then 
quit  the  iteration; 

else  if  ||r’|| 2  >  || r-* 1 1| 2  then  {  stage  switching} 

j  =  j  +  1  {  Update  the  stage  index} 

X  =  j  x  s  {  Update  the  size  of  finalist} 

else 

Fk  =  F  {  Update  the  finalist} 
rk  =  r  {  Update  the  residue} 
k  =  k  +  1 

end  if 

until  halting  condition  true; 

Output:  x p  =  T  }  y  {  Prediction  of  non-zero  coefficients} 
Algorithm  1:  Sparsity  adaptive  matching  pursuit  (SAMP) 


tions,  our  proposed  sparsity  adaptive  matching  pursuit  (SAMP)  is 
designed  to  take  advantages  of  both  bottom-up  and  top-down  ap¬ 
proaches. 

Fig.  3  shows  the  conceptual  diagram  of  the  SAMP  in  the  k- 
th  iteration.  One  can  observe  that  it  is  quite  similar  to  that  of  the 
SP/CoSaMP  algorithms  in  Fig.  1(c)  except  that  the  sizes  of  candi¬ 
date  set  |Cfc|  and  finalist  |Afc|  are  adaptive.  This  key  innovation  en¬ 
ables  the  SAMP  to  conduct  blind  recovery  without  priori  information 
of  K.  The  optimal  values  for  \Ck\  and  | Fk\  in  each  iteration  remain 
as  open  questions.  For  simplicity,  we  divide  the  recovery  process 
into  several  stages,  each  of  which  contains  several  iterations.  |Afc| 
is  kept  fixed  for  iterations  in  the  same  stage  and  increased  by  a  step 
size  s  <  K  between  two  consecutive  stages.  Also,  just  as  in  the  SP. 
the  candidate  set  is  chosen  as  \Ck\  =  2 1 Fk \ . 

Algorithm  1  presents  the  pseudo  code  of  the  SAMP.  Here, 
X  =  | Fk\  represents  the  size  of  finalist  and  for  a  vector  a,  function 
Max(a,X)  returns  X  indices  corresponding  to  the  largest  absolute 
values  of  a.  Also,  for  a  set  A  €  {1,  ■  •  •  ,  N},  <f>A  is  the  submatrix 
of  <E>  with  indices  indices  i  £  A.  At  the  fc-th  iteration,  Sk,  Ck  Fk, 
Vk  represent  the  short  list,  the  candidate  list,  the  finalist  and  the 
observation  residual,  respectively.  For  practical  applications,  two 
immediate  questions  about  Algorithm  1  are:  (1)  What  are  the  halting 
conditions?  (2)  How  to  choose  the  step  size  s? 

Halting  conditions :  Just  as  in  the  SP,  for  sparse  signals,  the 


SAMP  stops  when  the  residual’s  norm  || r || 2  is  smaller  than  a  cer¬ 
tain  threshold  e.  Here,  e  =  0  for  noiseless  measurements  and  e  can 
be  chosen  as  the  noise  energy  for  noisy  measurements.  Halting  con¬ 
dition  for  compressible  signals  is  more  complicated.  In  this  case, 
there  is  no  known  optimal  way  to  stop  the  algorithm,  even  with  con¬ 
vex  relaxation  algorithms.  One  common  approach  is  to  halt  when  a 
relative  residue  improvement  between  two  consecutive  iterations  is 
smaller  than  a  certain  threshold.  The  underlying  intuition  is  that  it 
would  not  worth  to  take  more  costly  iterations  if  the  resulting  im¬ 
provement  is  too  small.  For  example,  in  the  GPSR  algorithm  of  [3], 
the  program  stops  when  coordinates  in  the  finalist  changes  by  a  rela¬ 
tive  amount  less  than  a  threshold.  Based  on  this  principle,  we  suggest 
that  the  SAMP  halts  when  the  relative  change  of  reconstructed  sig¬ 
nal’s  energy  between  two  consecutive  stages  is  smaller  than  a  certain 
threshold. 

The  step  size  s:  The  SAMP  algorithm  only  requires  s  <  K.  To 
avoid  overestimation,  the  safest  choice  is  certainly  s  =  1  if  K  is 
unknown.  However,  there  is  a  trade-off  between  s  and  the  recovery 
speed  as  smaller  s  requires  more  iterations.  Also,  as  we  will  show 
in  Section  5,  the  choice  of  s  also  depends  on  the  magnitude  distri¬ 
bution  of  the  input  signal.  Empirical  results  suggest  that  small  s  is 
preferable  for  signal  with  (exponentially)  decayed  magnitude,  while 
large  s  is  advantageous  for  binary  sparse  signal.  The  derivation  of 
the  optimal  value  for  s  remains  as  an  open  question. 

3.2.  SAMP  vs.  existing  greedy  algorithms 

From  practical  perspective,  the  most  prominent  feature  of  the  SAMP 
lies  in  the  fact  that  it  does  not  require  K  as  an  input  parameter.  Front 
the  theoretical  point  of  view,  it  still  offers  a  strong  guarantee  for  exact 
recovery  and  stability,  as  we  will  show  in  Section  4.  Also,  just  as  the 
STOMP,  the  SAMP  adopts  a  stagewise  approach  to  expand  the  true 
support  set  stage  by  stage.  In  the  meantime,  it  takes  the  advantage 
of  the  backtracking  idea  in  the  SP/CoSAMP  to  refine  the  estimate  of 
true  support  set  at  each  iteration.  In  this  light,  it  is  a  combination  of 
both  bottom-up  and  top-down  principles. 

We  also  want  to  point  out  that  the  SAMP  provides  a  genarl- 
ized  framework  for  the  OMP  and  the  SP.  Note  that  when  s  =  1, 
SAMP  can  be  roughly  regarded  as  the  (generalized)  OMP  associ¬ 
ated  with  refinement  feature  that  can  remove  bad  coordinates  during 
iterations.  In  this  case,  the  SAMP  is  always  more  accurate  than  the 
OMP  although  it  may  require  a  few  more  iterations  to  achieve  that 
accuracy.  In  addition,  when  s  =  K ,  SAMP  becomes  exactly  SP 
if  the  restricted  isometry  property  (RIP)  condition  of  measurement 
matrix  is  satisfied.  In  this  case,  it  only  needs  one  stage  to  find  the  A'- 
sparse  approximation  of  the  signal.  Even  when  s  <  K,  each  stage 
in  the  SAMP  still  uses  a  similar  principle  of  the  SP,  i.e.  identifying  a 
portion  of  coordinates  in  the  true  support  set  and  then  using  several 
iterations  to  refine  this  estimate.  However,  in  general,  SAMP  and 
SP  behave  differently.  Compared  with  the  SP,  SAMP  establishes  a 
different  goal:  at  each  stage  it  attempts  to  discover  a  smaller  number 
of  coordinates  in  the  true  support  but  expects  a  higher  accuracy. 


4.  THEORETICAL  PERFORMANCE  ANALYSIS 

This  section  describes  our  theoretical  analysis  of  the  behavior  of 
SAMP  for  sparse  and  compressible  signals  in  both  noiseless  and 
noisy  cases.  Because  the  proofs  are  mainly  based  on  the  proof 
framework  of  SP,  the  following  theorems  are  formatted  in  parallel 
with  those  in  [7] 


Theorem  1  ( Exact  recovery  for  sparse  signals ):  Assume  x  £  R  v 
is  a  A'-sparse  signal  and  the  corresponding  measurement  y  =  <E>  x. 
Let  Ks  =  s\K/s\.  If  the  sensing  matrix  <1?  satisfies  the  RIP  with 
parameter:  <$3 ks  <  0.06,  SAMP  algorithm  guarantees  exact  recov¬ 
ery  of  x  from  y  via  a  finite  number  of  iterations 

Proof  Draff.  The  proof  is  mainly  based  on  the  following  lemma: 

Lemma  1 :  If  the  sensing  matrix  satisfies  the  above  condition  of  RIP: 

•  The  stage  [A/s]  is  equivalent  to  SP  algorithm  with  estimated 
sparsity  A's,  except  that  they  have  different  initial  final  lists 
and  initial  observation  data  vectors 

•  SAMP  recovers  the  target  signal  exactly  after  completing 
stage  [A/s] 

Proof  of  Lemma  1:  At  the  stage  [A/s],  both  finalist  and  short  list 
have  size  of  Ks  >  A.  Given  the  sizes  of  these  lists,  preliminary 
test  and  final  test  of  SAMP  are  similar  to  those  of  SP  with  the  cor¬ 
responding  value  of  sparsity.  The  only  difference  is  that  while  SP 
algorithm  has  empty  initial  finalist  and  full  initial  observation  data, 
stage  [A'/s]  has  the  finalist  and  observation  residual  of  the  last  iter¬ 
ation  as  its  initialization.  This  is  the  first  part  of  Lemma  1 . 

The  second  part  of  Lemma  1  is  derived  from  the  fact  that  the 
condition  of  convergence  of  SP  algorithm  in  [7]  does  not  depend 
on  the  those  initial  values  but  the  preliminary  test  and  final  test.  In 
particular,  it  is  based  on  the  following  observations  : 

•  Energy  of  the  part  of  signal  x  not  captured  by  the  current  fi¬ 
nalist  is  a  constant  times  smaller  than  that  of  signal  x  not 
captured  by  the  finalist  in  previous  iteration. 

•  Energy  of  observation  residual  of  the  current  iteration  is  a 
constant  times  smaller  than  that  of  previous  iteration 

When  the  condition  of  RIP  is  satisfied,  both  above  constants  are 
smaller  than  one  which  results  in  the  exact  recovery  after  a  finite 
number  of  iterations.  This  is  the  main  content  of  Theorem  2  and 
Theorem  7  in  [7],  Because  the  last  stage  is  equivalent  to  SP  with 
estimated  sparsity  level  A's,  it  is  obvious  that  the  target  signal  will  be 
exactly  recovered  after  this  stage  if  the  condition  of  RIP  of  parameter 
Ks  is  satisfied  .  To  complete  the  proof,  it  is  sufficient  to  show  that  the 
SAMP  algorithm  never  gets  stuck  at  any  iteration  of  any  stage,  i.e. 
it  takes  a  finite  number  of  iterations  up  to  the  stage  [A'/s].  Because 
at  each  stage,  the  finalist  (whose  size  is  assumed  to  be  p ),  add  and 
discard  some  coordinates  and  there  is  a  finite  number  of  coordinates, 
there  are  a  finite  number  of  combinations,  at  most  (^),  where  N  is 
the  length  of  the  signal.  Thus,  if  there  were  an  infinite  number  of 
iterations,  final  lists  would  be  repeated.  But  this  is  contradict  with 
the  stage  switching  condition  that  the  observation  residual  is  always 
monotonic  decreasing.  Hence,  Theorem  1  follows. 

Theorem  2  (Stability  for  sparse  signals ):  With  the  same  assumption 
and  notation  of  Theorem  1  and  assume  the  measurement  vector  is 
contaminated  with  noise:  y  =  <I>  x  +  e.  Let  energy  of  noise  be  o2. 
If  the  sensing  matrix  satisfies  the  RIP  with  parameter:  53ks  <  0.06, 
the  signal  approximation  x  of  SAMP  algorithm  satisfies: 

II*  —  *||  2  <  Cks<7  (2) 

where  cKs  =  ( 1  +  53ks ) / (<WS  (1  -  53ks  )) 

Theorem  3  (Stability  for  compressible  signals):  Assume  when  the 
algorithm  stops,  the  number  of  coordinates  in  the  finalist  is  Kstp. 
With  the  same  assumption  of  Theorem  2,  if  the  sensing  matrix  satis¬ 
fies  the  RIP  with  parameter:  5eKstp  <  0.03,  the  signal  approxima¬ 
tion  x  of  SAMP  algorithm  satisfies: 

||*  -  *||2  <  C2Kstp  (<7+^(1  +  S6K3tp)/Kstp\\x  -  XKstp  ||l)  (3) 


Similarly,  the  proof  of  Theorem  2  and  Theorem  3  is  based  on 
Lemma  1  and  the  independence  of  corresponding  theorems  of  SP 
algorithm  with  its  initialization.  We  omit  the  detail  proof  due  to 
space  limitation. 

The  above  theorems  are  sufficient  conditions  of  SAMP  for  exact 
recovery  and  stability.  They  are  slightly  more  restrictive  than  corre¬ 
sponding  results  of  SP  algorithms  because  the  true  sparsity  level  K 
is  always  smaller  than  or  equals  the  estimated  one  Ks.  This  may  be 
regarded  as  an  additional  cost  for  not  having  precise  information  of 
sparsity.  On  the  other  hand,  the  proofs  also  show  that  these  sufficient 
conditions  may  not  optimal  or  tight  enough  because  they  only  con¬ 
sider  the  last  stage  and  ignore  the  influence  of  previous  stages  to  the 
whole  performance.  This  issue  is  one  of  our  future  works. 

5.  SIMULATION  RESULTS 

This  section  compares  the  simulation  results  of  the  proposed  SAMP 
with  other  greedy  algorithms  and  the  h  optimization  algorithm.  We 
also  observe  some  interesting  performance  behaviors  that  could  not 
be  justified  by  our  theoretical  analysis,  especially  when  measure¬ 
ments  are  insufficient  for  exact  recovery.  These  results  imply  the 
limitations  of  the  sufficient  conditions  presented  in  Section  4.  Some 
heuristic  arguments  are  presented  to  complement  the  theoretical  part 
and  demystify  these  observation  results. 

5.1.  Experiment  1 

In  this  experiment,  the  signals  of  interests  are  Gaussian  or  binary 
sparse  signals  with  length  of  N  =  256.  The  partial  FFT  sensing 
operator  is  used  with  a  fixed  number  of  measurements  M  =  128. 
Our  aim  is  to  investigate  the  probability  of  exact  reconstruction  vs. 
the  signal  sparsity  A  for  a  given  M.  Different  sparsitys  are  chosen 
from  K  =  10  to  K  =  70  and  for  each  A',  500  simulations  were 
conducted  to  calculate  the  probabilities  of  exact  reconstruction  for 
different  algorithms.  Fig.  4(a)  and  Fig.  4(b)  demonstrate  the  results 
for  Gaussian  sparse  and  binary  sparse  signals,  respectively. 

As  can  be  seen,  for  Gaussian  sparse  signals,  performance  of 
the  SAMP  far  exceeds  that  of  all  other  algorithms,  including  the  l\- 
minimization.  While  all  other  algorithms  start  to  fail  when  sparisty 
K  >  40,  the  SAMP  still  can  afford  until  sparsity  K  >  60 — nearly 
equal  a  half  of  the  number  of  measurements.  However,  for  binary 
sparse  signals,  the  SAMP,  along  with  SP,  CoSaMP,  are  worse  than 
li  -minimization.  They  start  to  fail  when  sparisty  A  >  30  while 
li  -minimization  begins  to  fail  at  A  >  40 

5.2.  Experiment  2 

This  experiment  investigates  the  probability  of  exact  recovery  vs. 
the  number  of  measurements,  given  a  fixed  signal  sparsity  K.  We 
use  the  same  setups  of  experiment  1  and  choose  K  =  20,  M  £ 
(50,  60,  70,  80,  90, 100).  For  each  value  of  M,  we  generate  a  signal 
x  of  sparsity  K  and  its  measurements  y  =  4>x.  Then  we  use  above 
algorithms  to  recover  x.  This  procedure  is  repeated  500  times  for 
each  value  of  M.  We  then  calculate  the  probabilities  of  exact  recon¬ 
struction.  Fig.  5(a)  and  Fig.  5(b)  depict  these  probability  curves  of 
Gaussian  and  sparse  signals,  respectively.  The  numerical  values  on 
x-axis  denote  the  number  of  measurements  M  and  those  on  y-axis 
represent  probability  of  exact  recovery. 

Again,  we  see  that  SAMP  and  h -minimization  are  the  best  algo¬ 
rithms  for  recovering  Gaussian  and  sparse  signals,  respectively.  It  is 
also  interesting  to  observe  that  when  the  number  of  measurements  is 
insufficient  for  guarantee  of  exact  recovery,  the  probability  of  exact 
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Fig.  4.  Prob.  of  exact  recovery  vs.  the  signal  sparsity  K.  Here,  the 
test  signal  is  of  length  N  =  256  and  the  number  of  measurements 
is  fixed  as  M  =  128.  (a)  Gaussian  sparse  signal,  (b)  Binary  sparse 
signal. 


recovery  of  SAMP  depends  on  its  step  size  and  signal  types.  In  par¬ 
ticular,  for  Gaussian  sparse  signals,  SAMP  with  a  smaller  step  size 
gets  a  higher  chance  of  recovering  signals  exactly,  given  the  same 
number  of  measurements.  On  the  contrary,  for  binary  sparse  signals, 
SAMP  with  a  bigger  step  size  (e.g  SP)  gets  a  better  chance  of  exact 
recovery  of  signals.  Although  these  observation  could  not  be  justi¬ 
fied  by  theorems  of  sufficient  conditions,  they  may  be  heuristically 
justified  as  follows. 

When  a  signal  is  exponentially  decayed,  a  preliminary  test  which 
is  based  on  the  principle  of  maximal  correlation,  is  not  accurate  if  a 
large  number  of  coordinates  are  admitted  into  the  short  list.  It  means 
that  many  of  them  might  be  the  wrong  coordinates(not  in  the  true 
support  set).  Thus,  for  this  type  of  signal,  select  few  but  with  more 
caution  at  each  stage  is  more  efficient  and  accurate.  As  a  result, 
SAMP  with  a  smaller  step  size  is  more  accurate  than  SP. 

On  the  contrary,  when  a  signal  is  binary  sparse,  a  preliminary 
test  is  still  decently  accurate  even  if  we  admit  more  coordinates  into 
the  short  list.  However,  to  justify  why  SAMP  with  a  smaller  step  size 
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Fig.  5.  Prob.  of  exact  recovery  vs.  the  number  of  measurements 
M.  Here,  the  test  signal  is  of  length  N  =  256  and  the  number  of 
measurements  is  fixed  as  M  =  128.  (a)  Gaussian  sparse  signal,  (b) 
Binary  sparse  signal. 


works  worse  in  this  case,  it  is  necessary  to  look  into  the  structure  of 
its  final  test.  At  the  stage  k,  we  first  find  the  least  square  signal  ap¬ 
proximation  in  the  subspace  spanned  by  columns  whose  coordinates 
are  in  the  candidate  set  Ck ■  Then  we  admit  a  subset  of  those  co¬ 
ordinates  whose  entries  are  largest  in  magnitude  to  be  the  final  list. 
This  process  expects  to  capture  the  coordinates  in  T  PI  Ck,  where 
T  denotes  the  true  support  set.  If  T  C  Ck,  the  process  is  likely  to 
obtain  it  expectation.  However,  when  T  \  Ck  is  not  empty,  it  results 
in  incoherent  noise  yr\ck  that  in  turn,  distorts  the  least  square  sig¬ 
nal  approximation.  Our  hypothesis  is  that  energy  of  this  incoherent 
noise  \\yT\ck  II2  would  affect  the  accuracy  of  the  final  list  Fk  specif¬ 
ically  and  degrade  the  whole  performance  in  general.  SAMP  with 
smaller  step  size  s  is  less  efficient  because  size  of  candidate  list  Ck 
which  is  proportional  to  s  is  small.  Thus,  even  Ck  C  T,  ||t/T\cfc  II 2  is 
still  large  and  that  results  in  very  high  incoherent  noise.  On  the  other 
hands,  SP  is  more  efficient  because  its  candidate  list  \Ck\  =  A'  is 
large  and  due  to  the  efficiency  of  the  preliminary  test,  \\yTnck  H2  is 
also  large  and  thus, the  incoherent  noise  energy  ||t/T\cJ|2  is  rela- 


tively  small 

Finally,  Fig.  6  depicts  the  stagewise  recovery  and  its  incoherent 
noise  VT\ck  =  &T\ckxT\ck  for  binary  sparse  and  decayed  sparse 
signals,  respectively.  Due  to  RIP  of  <f>T\cfc>  II  J/r\c?fc  II 2  is  propor¬ 
tional  to  ||a;T\cfc  H2.  These  figures  demonstrate  that  energy  of  inco¬ 
herent  noise  when  signal  is  binary  sparse  is  larger  than  when  signal 
is  rapidly  decayed.  In  other  words,  stagewise  recovery  of  SAMP 
is  more  efficient  with  rapidly  decayed  signal  and  less  efficient  with 
binary  sparse  signal. 
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Fig.  6.  Incoherent  noise  generated  by  stagewise  recovery,  (a)  Binary 
sparse  signal,  (b)  Decayed  sparse  signal. 


5.3.  Experiment  3 

In  this  experiment  set,  we  compare  the  performance  of  the  SAMP, 
the  StOMP.  the  SP  and  the  h -minimization  in  the  practical  large- 
scale  compressed  imaging  scenario.  The  software  GPSR  is  used  for 
h -minimization  because  of  its  fast  speed  and  good  performance  [3], 
In  addition,  for  the  SP  algorithm,  the  input  parameter  K  is  estimated 
by  setting  it  equal  to  the  number  of  nonzero  entries  that  the  SAMP 
can  detect. 

Three  512  x  512  test  images  Lena ,  Barbara  and  Boat  were  cho¬ 
sen  and  the  sparsifying  matrix  is  the  popular  Daubechies  wavelet 
9/7  wavelet.  Structurally  random  matrices  were  used  as  the  sam¬ 
pling  operator  due  to  their  fast  and  efficient  implementations  [9]. 
Table.  1  summarizes  the  PSNR  results  of  different  algorithms  and 
Fig.  7  shows  the  visual  reconstructions  of  the  image  Lena  512  x  512 
from  M  =  N/ 5  measurements.  The  above  table  and  figures  imply 
that  in  the  practical  compressed  sensing,  performance  of  the  pro¬ 
posed  SAMP  is  comparable  to  that  of  Linear  Programming  and  ex¬ 
ceeds  several  other  greedy  algorithms. 

In  general,  the  SAMP  may  require  more  iterations  than  other 
greedy  algorithms  such  as  StOMP,  SP  and  CoSaMP,  especially  when 
the  step  size  s  is  much  smaller  than  the  true  sparisty  K .  In  the  ex¬ 
treme  case,  when  step  size  s  "as  1,  SAMP  becomes  the  generalized 
OMP  and  it  would  require  at  least  K  iterations.  As  depicted  in  the 
experiment  1  and  experiment  2,  with  compressible  sparse  signals, 
when  we  increase  step  size  s,  SAMP  takes  fewer  iterations  but  its 


performance  is  gradually  degraded.  This  is  a  trade-off  between  com¬ 
putational  complexity  and  performance.  How  to  define  optimal  val¬ 
ues  of  step  size  s,  given  some  prior  model  of  compressible  signals  is 
our  future  research  question. 

6.  CONCLUSIONS 

In  this  paper,  a  novel  greedy  pursuit  algorithm,  called  the  sparsity 
adaptive  Matching  Pursuit,  is  proposed  and  analyzed  for  reconstruc¬ 
tion  applications  in  compressed  sensing.  As  its  name  suggests,  this 
reconstruction  algorithm  is  most  featured  of  not  requiring  informa¬ 
tion  of  sparsity  of  target  signals  as  a  prior.  It  not  only  releases  a  com¬ 
mon  limitation  of  existing  greedy  pursuit  algorithms  but  also  keeps 
performance  comparable  with  that  of  strongest  algorithms  such  as 
SP,  CoSaMP  or  linear  programming.  The  underlying  intuition  of 
SAMP  which  is  similar  to  that  of  the  EM  algorithm  is  to  alterna¬ 
tively  estimate  the  sparsity  and  the  true  support  set.  Extensive  ex¬ 
periment  results  confirm  that  SAMP  is  very  appropriate  for  recon¬ 
structing  compressible  sparse  signal  where  its  magnitudes  are  de¬ 
cayed  rapidly. 
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Table  1.  Comparison  of  algorithms’  objective  performance:  PSNR  in  dB 


M/N 

(Sampling  rate) 

Lena 

Boat 

Barbara 

GPSR 

StOMP 

SP 

SAMP 

GPSR 

StOMP 

SP 

SAMP 

GPSR 

StOMP 

SP 

SAMP 

0.1 

25.90 

22.63 

23.43 

23.52 

18.40 

20.77 

20.91 

0.2 

28.48 

25.73 

23.58 

25.87 

25.86 

22.58 

21.24 

22.63 

22.80 

0.3 

32.07 

26.88 

28.68 

28.85 

25.04 

23.07 

23.58 

24.87 

0.4 

33.99 

31.06 

27.45 

26.25 

25.83 

27.80 

0.5 

35.78 

35.59 

34.73 

35.42 

32.37 

33.33 

32.71 

30.11 

30.27 

29.33 

30.19 

Fig.  7.  Reconstructed  512  x  512  Lena  images  front  M/N  =  20%  sampling  rate.  Results  of  (a)  GPSR:  28.65dB;  (b)  StOMP:  26.29dB;  (c) 
SP:  28.47dB;  (d)SAMP:  28.48dB 


